Section 4 Scenario Analysis: WCS Priority Landscapes

# turn off the s2 processing 
## https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
sf::sf_use_s2(FALSE)

See the USFS Wildfire Crisis Strategy (WCS) document for background information.

4.1 Scenarios Considered

There are three scenarios considered in this analysis of the constraints to mechanical forest health and fuel reduction treatments:

# mechanical mgmt allowed if:
n_sc_temp <- 3
data.frame(
  scenario = 1:3
  , cover = rep("Forest or Shrubland",n_sc_temp)
  , protected = rep("Not Protected",n_sc_temp)
  , slope = c("<40%","<60%","<60%")
  , road = c("<1,000ft","<2,000ft","<2,000ft")
  , riparian = c(">100ft",">100ft",">50ft")
  , admin = c("No Designation","No Designation","Any Designation")
) %>% 
  tidyr::pivot_longer(
    cols = -c(scenario)
  ) %>% 
  tidyr::pivot_wider(
    names_from = scenario
    , values_from = value
    , names_prefix = "scenario_"
  ) %>% 
  dplyr::mutate(
    name = factor(
      name
      , levels = c(
        "cover"
        , "protected"
        , "slope"
        , "road"
        , "riparian"
        , "admin"
      )
      , labels = c(
        "NLCD Cover Type"
        , "Protected or<br>IRA Status"
        , "Slope"
        , "Distance to<br>Nearest Road"
        , "Riparian<br>Buffer"
        , "Administrative<br>Designation"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::arrange(name) %>% 
  dplyr::mutate(
    lvl = paste0("L",dplyr::row_number()-1,":")
  ) %>% 
  dplyr::relocate(lvl) %>% 
# make table
kableExtra::kable(
    caption = "Scenarios Considered for Determining Mechanical Treatment Operability<br>Mechanical treatment allowed if:"
    , col.names = c(
      ""
      , ""
      , "Scenario 1<br>Most Constrained"
      , "Scenario 2<br>Moderate"
      , "Scenario 3<br>Least Constrained"
    )
    , escape = F
  ) %>%  
  kable_classic(full_width=T) %>% 
  kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE)

Table 4.1: Table 4.2: Scenarios Considered for Determining Mechanical Treatment Operability
Mechanical treatment allowed if:

		Scenario 1<br>Most Constrained	Scenario 2<br>Moderate	Scenario 3<br>Least Constrained
L0:	NLCD Cover Type	Forest or Shrubland	Forest or Shrubland	Forest or Shrubland
L1:	Protected or<br>IRA Status	Not Protected	Not Protected	Not Protected
L2:	Slope	<40%	<60%	<60%
L3:	Distance to<br>Nearest Road	<1,000ft	<2,000ft	<2,000ft
L4:	Riparian<br>Buffer	>100ft	>100ft	>50ft
L5:	Administrative<br>Designation	No Designation	No Designation	Any Designation

The figure below illustrates the workflow used to quantify the amount of land available for mechanical forest health and risk reduction fuel treatments by considering layered operational constraints. This example illustrates scenario 2 for a 10,252 hectare (25,000 acre) fireshed project area near Rustic, Colorado in the Colorado Front Range priority landscape.

4.2 Data Preparation

4.2.1 Load WCS landscapes

WCS priority landscape spatial data from the USFS Geospatial Data Discovery site (accessed 2023-09-13).

# load data
wf_landscapes <- sf::read_sf("../data/Wildfire_Crisis_Strategy_Landscapes/Wildfire_Crisis_Strategy_Landscapes_(Feature_Layer).shp") %>% 
  dplyr::rename_with(tolower) %>% 
  sf::st_make_valid() %>% 
  dplyr::mutate(
    hectares = (as.numeric(sf::st_area(geometry))/10000) 
    , Mil.Hectares = hectares %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
    , acres = (as.numeric(sf::st_area(geometry))/4046.85642) 
    , Mil.Acres = acres %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
  ) %>% 
  dplyr::left_join(
    data.frame(state = datasets::state.name, state_abb = datasets::state.abb)
    , by = join_by("state")
  ) %>% 
  dplyr::mutate(
    area_name = paste0(
      ifelse(is.na(state_abb) | state_abb == "", "UNK", state_abb)
      , ": "
      , name
    )
  ) %>% 
  # manual label placement
  dplyr::left_join(
    readr::read_csv("../data/wildfirepriority_WCS202308/area_label_placement.csv")
    , by = dplyr::join_by("area_name")
  )
#rename sf geom column
  names(wf_landscapes)[names(wf_landscapes)==tolower(attr(wf_landscapes, "sf_column"))] = "geometry"
  sf::st_geometry(wf_landscapes) = "geometry"
# set crs
  transform_crs <- sf::st_crs(wf_landscapes)
# view
  wf_landscapes %>% dplyr::glimpse()

4.2.2 Load landscape-level constraint data

Landscape-level constraint analysis data (i.e., row unique by WCS landscape) was created via this Google Earth Engine script.

# load landscape constraint scenario data
constrained_by_scnro_ls <-
  list.files("../data/wildfirepriority_WCS202308/landscapes/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_all_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority_WCS202308/landscapes/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::group_by(area_name,scenario_id) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # scenario description
    , scenario_lab = factor(
      scenario_id
      , levels = 1:3
      , labels = paste0("Scenario ", 1:3)
      , ordered = T
    ) %>% forcats::fct_rev()
    , scenario_desc = factor(
      scenario_id
      , levels = 1:3
      , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
      , ordered = T
    )
  )

4.2.3 Load fireshed project area data

The fireshed registry spatial data was obtained from the USFS Geospatial Data Discovery tool for firesheds and project areas (accessed 2023-05-03).

See this section for exploration of fireshed and project area spatial data

### fireshed spatial data
fireshed <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Fireshed/Fireshed_Registry%3A_Fireshed_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "area_ha"
      , "fireshed_id"
      , "fireshed_name"
      , "fireshed_code"
      , "fireshed_state"
      , "nopas"
      , "objectid"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "crisis_strategy"
      , "key_preformance_indicator"
      , "national_usfs_rank"
      , "national_all_land_rank"
      , "regional_usfs_rank"
      , "regional_all_land_rank"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
    # there is also a national_all_land_rank column
    , ntllandrank_pct_rank = dplyr::percent_rank(-national_all_land_rank)
    , ntllandrank_pct_rank_grp = dplyr::case_when(
        ntllandrank_pct_rank >= 1-0.01 ~ "Top 1%"
        , ntllandrank_pct_rank >= 1-0.05 ~ "Top 5%"
        , ntllandrank_pct_rank >= 1-0.10 ~ "Top 10%"
        , ntllandrank_pct_rank >= 1-0.25 ~ "Top 25%"
        , TRUE ~ "Bottom 75%"
      ) %>% 
      factor(
        levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
        , ordered = T
      )
    , crisis_strategy = ifelse(is.na(crisis_strategy),"Not High Risk",crisis_strategy) %>% 
      as.factor() %>% 
      forcats::fct_shift()
  )
  #rename sf geom column
    names(fireshed)[names(fireshed)==tolower(attr(fireshed, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed) = "geometry"
    # calculate area
    fireshed <- fireshed %>% 
      dplyr::mutate(
        fireshed_area_ha = as.numeric(sf::st_area(geometry))/10000
        , fireshed_area_acres = (fireshed_area_ha*10000)/4046.85642
      )
## fireshed_proj_area spatial data
fireshed_proj_area <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Project_Area/Fireshed_Registry%3A_Project_Area_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "fireshed_id"
      , "pa_id"
      , "pa_area_ha"
      , "objectid"
      , "pa_id2"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "pctrecentlydisturbed"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
  )
  #rename sf geom column
    names(fireshed_proj_area)[names(fireshed_proj_area)==tolower(attr(fireshed_proj_area, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed_proj_area) = "geometry"
    # calculate area
    fireshed_proj_area <- fireshed_proj_area %>% 
      dplyr::mutate(
        pa_area_ha = as.numeric(sf::st_area(geometry))/10000
        , pa_area_acres = (pa_area_ha*10000)/4046.85642
      ) %>% 
      # JOIN WITH FIRESHED DATA
      dplyr::inner_join(
        fireshed %>%
          sf::st_drop_geometry() %>%
          dplyr::select(fireshed_id, crisis_strategy, exp_total
                        , exposure_pct_rank, exposure_pct_rank_grp
          ) %>% 
          dplyr::rename(exposure_total=exp_total) %>% 
          dplyr::rename_with(
            ~ paste0("fireshed_",.x)
            , -c(fireshed_id)
          )
        , by = dplyr::join_by(fireshed_id)
      ) %>%
      dplyr::select(pa_id,pa_area_ha
                    ,exp_total,exposure_pct_rank,exposure_pct_rank_grp
                    , tidyselect::starts_with("fireshed_")
      ) %>% 
      dplyr::rename(exposure_total=exp_total) %>%
      dplyr::rename_with(
        ~ paste0("pa_", .x, recycle0 = TRUE)
        , tidyselect::starts_with("exp")
      )

Fireshed project area constraint analysis data (i.e., row unique by WCS landscape + FPA) was created via this Google Earth Engine script

# load fireshed constraint scenario data
constrained_by_scnro_ls_pa <-
  list.files("../data/wildfirepriority_WCS202308/fireshed/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_fireshed_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority_WCS202308/fireshed/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,pa_id,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  dplyr::group_by(area_name,pa_id,scenario_id) %>%
  dplyr::arrange(desc(pct_pa_intrsct)) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::filter(pct_pa_intrsct>=0.00001) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # calculate level of constraint
      , cnstrnt_lvl = dplyr::case_when(
          -pct_rdctn_total > 0.8 ~ 1
          , -pct_rdctn_total >= 0.6 ~ 2
          , -pct_rdctn_total >= 0.0 ~ 3
        )
      , cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("high constraint", "med. constraint", "low constraint")
          , ordered = T
        ) %>% forcats::fct_rev()
      , rmn_cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("0–19% treatable", "20–40% treatable", ">40% treatable")
          , ordered = T
        ) %>% forcats::fct_rev()
      # scenario description
      , scenario_lab = factor(
          scenario_id
          , levels = 1:3
          , labels = paste0("Scenario ", 1:3)
          , ordered = T
        ) %>% forcats::fct_rev()
      , scenario_desc = factor(
        scenario_id
        , levels = 1:3
        , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
        , ordered = T
      )
  ) %>% 
  # join with fireshed data
  dplyr::inner_join(
    fireshed_proj_area %>%
      sf::st_drop_geometry() %>%
      dplyr::select(pa_id, tidyselect::starts_with("pa_exposure"), tidyselect::starts_with("fireshed_")) %>%
      dplyr::mutate(pa_id=as.character(pa_id))
    , by = dplyr::join_by(pa_id)
  )

4.3 Geographic Delineation Description

Timeline of the research utilized to motivate and support the creation of the WCS:

• 2018 August - Toward share stewardship across the landscapes: an outcome-based investment strategy created a hotspot map of the source of wildfire risk to communities showing that 80% of community exposure from NFS lands originates on 20% of the NFS area.
• 2019 February - Evers et al. analyzed community wildfire exposure from national forests by associating conditions that affect exposure in the areas where wildfires ignite to conditions where exposure likely occurs
• 2019 May - Ager et al. mapped community firesheds by creating a continuous smoothed surface of predicted structure exposure from all FSim ignitions (see also Palaiologou et al. 2019 and Ager et al. 2019)
• 2020 - Evers et al. release Fireshed Registry as a geospatial dashboard and decision tool
• 2021 May - Ager et al. 2021 official release of the Fireshed Registry as a geospatial information system that organizes landscape risk to developed areas into a planning framework (see figure below)
• 2022 January - US Department of Agriculture, Forest Service launches a 10-year Wildfire Crisis Strategy
• 2022 April - 10 initial landscapes selected for WCS prioritization
• 2023 January - additional 11 landscapes selected for WCS prioritization

Nested spatial framework for firesheds. Each scale has specific functionality in terms of the planning processes. Firesheds are the broad scale unit of prioritization, but project areas within them are also prioritized as part of implementation of treatments. Project areas are roughly the size that national forests use for conducting vegetation and fuel management projects. The relative variation among firesheds compared to variation within them controls the relative emphasis on selecting firesheds versus individual planning areas. Ager et al. 2021

#define basemap
states_temp <- USAboundaries::us_states(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs)
cities_temp <- USAboundaries::us_cities(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs) %>% 
  dplyr::filter(
    toupper(city) %in% c(
      "DENVER"
      , "TUSCON"
      , "SALT LAKE CITY"
      , "LAS VEGAS"
      , "SAN DIEGO"
      , "LOS ANGELES"
      , "FRESNO"
      , "SAN FRANCISCO"
      , "SACRAMENTO"
      , "PORTLAND"
      , "SEATTLE"
      , "ALBUQUERQUE"
      , "RENO"
      , "BOISE"
      , "HELENA"
      , "BILLINGS"
    )
    | (toupper(city) == "PHOENIX" & state_abbr == "AZ")
  )
# plot basemap
plt_basemap <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  geom_sf_text(
    data = states_temp
    , mapping = aes(label = stusps)
    , color = "gray66"
    , size = 3
  ) +
  geom_sf(
    data = cities_temp
    , shape = 1
    , size = 1
    , color = "gray44"
  ) +
  geom_sf_text(
    data = cities_temp
    , mapping = aes(label = city)
    , color = "gray44"
    , size = 2
    # , hjust = -0.15
    , hjust = 0.5
    , vjust = -0.25
  ) +
  coord_sf(expand = F) +
  theme_void()

# plot basemap simple
plt_basemap_simple <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  coord_sf(expand = F) +
  theme_void()

4.3.1 Fireshed risk to communities (2021)

ggplot() + 
  geom_sf(
    data = fireshed
    , mapping = aes(fill=ntllandrank_pct_rank_grp)
  ) + 
  geom_sf(
    data = USAboundaries::us_states() %>% 
      dplyr::filter(!stusps %in% c("AK","HI","PR")) %>% 
      sf::st_transform(transform_crs)
    , fill = NA
    , color = "gray88"
  ) +
  scale_fill_manual(values=c("firebrick","orange3","khaki","seagreen3","royalblue")) +
  coord_sf(expand = F) +
  labs(
    fill = "Fireshed exposure ranks"
  ) +
  theme_void() + 
  theme(
    legend.position = "top"
    , legend.direction = "horizontal"
    , legend.title = element_text(size = 8)
    , legend.text = element_text(size = 7)
  )

National map of the 7,688 firesheds created from community wildfire transmission data (Evers et al. 2020). The fireshed boundaries were created with a process that delineates hotspots of fire transmission to buildings in adjacent or nearby communities. See the Methods section for details on delineating firesheds. Figure 2 from Ager et al. 2021 (p. 7)

4.3.2 Wilfire Crisis Strategy 21 Priority Landscapes (2022-2023)

plt_fshed <- 
plt_basemap + 
  geom_sf(
    data = fireshed %>%
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(
      color = paste0(investment, " investment")
      # , linetype = paste0(investment, " investment")
    )
    , fill=NA
    # , color="black" 
    , lwd=0.7
  ) +
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray11", "navy")) + # midnightblue
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
    # , linetype = "Landscapes"
  ) +
  # ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )
# plot
plt_fshed

ggplot2::ggsave(
    filename = paste0("../data/plt_map_ls_fshed.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 8.5
    , units = "in"
    , dpi = "print"
    , bg = "white"
  )

The boundaries of the 21 priority landscapes roughly follow the boundaries of “firesheds” prioritized to reduce wildfire transmission to developed areas (Figure 1). Firesheds are geographic delineations with an average area of approximately 250,000 acres (~100,000 hectares) that were created to organize the landscape into units for managing wildfire risk to communities.

4.3.3 Wilfire Crisis Strategy Treatment Objective

plt_basemap + 
  geom_sf(
    data = fireshed %>%
      dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.7
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "navy")) +
  # ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

4.3.4 Fireshed Project Areas

name_temp <- "Enchanted Circle"
# pa map
plt_fshed_pa <-
ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>% 
          dplyr::select(pa_id)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = fireshed_crisis_strategy, color = "Project Area\nboundary")
    , alpha = 0.8
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.79) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# pa map w/in
plt_fshed_pa_incl <- ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>%
          dplyr::mutate(
            # is_25in = ifelse(pct_pa_intrsct>=0.25,"≥25%","<25%")
            is_25in = ifelse(pct_pa_intrsct>=0.25,"included","excluded")
          ) %>% 
          dplyr::select(
            pa_id
            , is_25in
          )
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = is_25in, color = "Project Area\nboundary")
    , alpha = 0.9
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = viridis::viridis(n=4, direction = -1, alpha = 0.9)[c(1,3)]) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Project Area\nInclusions"
    # fill = "Included in this analysis\nif ≥25% area w/in landscape"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.825) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# combine
cowplot::plot_grid(
    plotlist = list(
      NULL
      , plt_fshed_pa + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
      , plt_fshed_pa_incl + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
    )
    , nrow = 1
    , rel_widths = c(0.01,1,0.01,1,0.01)
    , labels = c("","A","","B","")
    , label_size = 17
    , label_fontface = "bold"
    , label_colour = "gray11"
    # , label_y = 0.99
  ) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_ex_spatial_strctr.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 5.8
    , units = "in"
    , dpi = "print"
  )

1. Nested within firesheds are project areas of approximately 10,000 ha in size which represent the geographic unit at which vegetation and fuel management projects are planned (Ager et al. 2021).

2. Fireshed project areas that did not have at least 25% of area within the boundary of the Strategy priority landscape were excluded from this analysis. This cutoff was implemented under the assumption that with less than 25% of area available, treatment alone would not substantially affect wildfire behavior across the fireshed (North et al. 2015). For each fireshed project area included in this analysis, the entire area within the project area boundary was used to calculate the percent of mechanically available acreage including area outside of the priority landscape boundary.

4.4 Landscape-level analysis

# filter fireshed project area
constrained_by_scnro_ls_pa <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(pct_pa_intrsct>=0.25)

4.4.1 Reduction Treatable Area Table

# scnum_temp<-2
  tbl_temp <- constrained_by_scnro_ls %>% 
    # dplyr::filter(scenario_id==scnum_temp) %>% 
    dplyr::mutate(
      dplyr::across(
        tidyselect::ends_with("_ha")
        , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
      )
      , dplyr::across(
        tidyselect::starts_with("pct_")
        , ~ scales::percent(.x, accuracy = .1)
      )
    ) %>% 
    dplyr::select(
      area_name
      , scenario_lab
      , covertype_area_ha
      , pct_covertype_area
      , pct_rdctn1_protected
      , pct_rdctn2_slope
      , pct_rdctn3_roads
      , pct_rdctn4_riparian
      , pct_rdctn5_administrative
      , rmn5_administrative_area_ha
      , pct_rmn5_administrative
    ) %>% 
    dplyr::arrange(area_name,desc(scenario_lab))
  # make table 
  kableExtra::kable(
      tbl_temp %>% dplyr::select(-c(area_name))
      # , caption = paste0("<b><font color=navy>Scenario "
      #                    ,scnum_temp
      #                    ,"</font></b><br>percent reduction of different types of constraints on mechanical treatment"
      #         )
      , caption = "<b>Wildfire Crisis Strategy 21 Priority Landscapes</b><br>percent reduction of different types of constraints on mechanical treatment"
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n(%)"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Mechanically\nAvailable (ha)"
        , "Mechanically\nAvailable (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Area of\nPriority Landscape"=2, "Constraint\nLeast Flexible to Most Flexible" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tbl_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")

Table 4.3: Table 4.4: Wildfire Crisis Strategy 21 Priority Landscapes
percent reduction of different types of constraints on mechanical treatment

	Area of Priority Landscape		Constraint Least Flexible to Most Flexible
	Forest+Shrub (ha)	|Forest+Shrub (%	|Protecte	|Slope Steepne	s |Road Dista	ce |Riparian Bu	fer |Administra	ive |Mechanically Available	(ha) |Mechanically Availab
AZ: 4FRI
Scenario 1	2,128k	90.4%	-12.4%	-7.5%	-26.7%	-10.2%	-15.5%	589k	27.7%
Scenario 2	2,128k	90.4%	-12.4%	-2.1%	-13.3%	-14.0%	-19.9%	813k	38.2%
Scenario 3	2,128k	90.4%	-12.4%	-2.1%	-13.3%	-9.9%	0.0%	1,325k	62.3%
AZ: Prescott
Scenario 1	275k	87.2%	-12.9%	-16.1%	-25.1%	-11.9%	-4.2%	82k	29.8%
Scenario 2	275k	87.2%	-12.9%	-3.3%	-13.2%	-18.6%	-6.8%	124k	45.2%
Scenario 3	275k	87.2%	-12.9%	-3.3%	-13.2%	-13.2%	0.0%	158k	57.4%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1	1,119k	91.7%	-9.3%	-19.7%	-49.8%	-5.8%	-1.1%	161k	14.4%
Scenario 2	1,119k	91.7%	-9.3%	-5.0%	-46.6%	-10.1%	-2.2%	300k	26.8%
Scenario 3	1,119k	91.7%	-9.3%	-5.0%	-46.6%	-7.3%	0.0%	357k	31.9%
CA: North Yuba
Scenario 1	139k	96.5%	-11.6%	-33.4%	-7.8%	-14.7%	-2.3%	42k	30.2%
Scenario 2	139k	96.5%	-11.6%	-13.9%	-3.2%	-23.0%	-3.3%	63k	45.0%
Scenario 3	139k	96.5%	-11.6%	-13.9%	-3.2%	-16.7%	0.0%	76k	54.6%
CA: Plumas Community Protection
Scenario 1	101k	93.3%	0.0%	-23.2%	-19.8%	-21.5%	-1.9%	34k	33.6%
Scenario 2	101k	93.3%	0.0%	-5.2%	-8.1%	-33.9%	-2.9%	50k	49.9%
Scenario 3	101k	93.3%	0.0%	-5.2%	-8.1%	-24.4%	0.0%	63k	62.3%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1	1,408k	86.0%	-57.0%	-21.0%	-8.3%	-2.3%	-1.4%	139k	9.9%
Scenario 2	1,408k	86.0%	-57.0%	-9.6%	-6.5%	-4.2%	-3.0%	278k	19.7%
Scenario 3	1,408k	86.0%	-57.0%	-9.6%	-6.5%	-3.0%	0.0%	336k	23.9%
CA: Stanislaus
Scenario 1	116k	94.3%	-2.9%	-26.7%	-10.6%	-24.7%	-0.6%	40k	34.5%
Scenario 2	116k	94.3%	-2.9%	-10.7%	-3.9%	-35.8%	-0.8%	53k	46.0%
Scenario 3	116k	94.3%	-2.9%	-10.7%	-3.9%	-26.1%	0.0%	66k	56.5%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1	554k	83.5%	-21.4%	-37.8%	-10.4%	-8.7%	-9.2%	69k	12.4%
Scenario 2	554k	83.5%	-21.4%	-12.4%	-7.8%	-17.3%	-18.5%	126k	22.7%
Scenario 3	554k	83.5%	-21.4%	-12.4%	-7.8%	-12.5%	0.0%	254k	45.9%
CO: Colorado Front Range
Scenario 1	1,052k	72.6%	-20.4%	-19.9%	-20.0%	-11.3%	-1.6%	281k	26.7%
Scenario 2	1,052k	72.6%	-20.4%	-5.6%	-12.0%	-17.3%	-3.2%	436k	41.4%
Scenario 3	1,052k	72.6%	-20.4%	-5.6%	-12.0%	-12.3%	0.0%	522k	49.6%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1	700k	89.0%	-43.0%	-19.7%	-9.1%	-1.4%	-0.3%	185k	26.5%
Scenario 2	700k	89.0%	-43.0%	-6.7%	-5.5%	-2.3%	-0.7%	293k	41.8%
Scenario 3	700k	89.0%	-43.0%	-6.7%	-5.5%	-1.6%	0.0%	302k	43.2%
ID: Southwest Idaho
Scenario 1	611k	87.6%	-13.3%	-28.2%	-10.4%	-2.0%	-0.3%	279k	45.8%
Scenario 2	611k	87.6%	-13.3%	-7.0%	-6.6%	-2.8%	-0.8%	425k	69.5%
Scenario 3	611k	87.6%	-13.3%	-7.0%	-6.6%	-2.0%	0.0%	435k	71.2%
MT: Kootenai Complex
Scenario 1	372k	91.4%	-21.4%	-15.0%	-7.8%	-6.9%	-10.6%	142k	38.3%
Scenario 2	372k	91.4%	-21.4%	-3.6%	-2.1%	-8.6%	-13.8%	188k	50.6%
Scenario 3	372k	91.4%	-21.4%	-3.6%	-2.1%	-6.1%	0.0%	249k	66.9%
NM: Enchanted Circle
Scenario 1	506k	85.6%	-8.5%	-21.5%	-33.4%	-7.5%	-2.2%	136k	26.8%
Scenario 2	506k	85.6%	-8.5%	-5.2%	-25.7%	-11.7%	-3.7%	229k	45.2%
Scenario 3	506k	85.6%	-8.5%	-5.2%	-25.7%	-8.2%	0.0%	265k	52.3%
NV: Sierra and Elko Fronts
Scenario 1	1,001k	73.1%	-26.6%	-9.8%	-26.8%	-6.2%	-0.9%	297k	29.6%
Scenario 2	1,001k	73.1%	-26.6%	-2.0%	-15.7%	-8.9%	-1.5%	453k	45.3%
Scenario 3	1,001k	73.1%	-26.6%	-2.0%	-15.7%	-6.3%	0.0%	495k	49.4%
OR: Central Oregon
Scenario 1	959k	87.1%	-7.0%	-2.8%	-15.3%	-5.4%	-10.0%	571k	59.5%
Scenario 2	959k	87.1%	-7.0%	-0.6%	-5.0%	-6.4%	-11.3%	669k	69.7%
Scenario 3	959k	87.1%	-7.0%	-0.6%	-5.0%	-4.5%	0.0%	796k	83.0%
OR: Klamath River Basin
Scenario 1	2,675k	77.0%	-19.3%	-14.9%	-17.7%	-5.4%	-4.4%	1,023k	38.3%
Scenario 2	2,675k	77.0%	-19.3%	-4.9%	-7.9%	-8.3%	-7.8%	1,386k	51.8%
Scenario 3	2,675k	77.0%	-19.3%	-4.9%	-7.9%	-5.9%	0.0%	1,659k	62.0%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1	326k	76.3%	-20.1%	-14.3%	-17.1%	-9.0%	-15.7%	78k	23.8%
Scenario 2	326k	76.3%	-20.1%	-4.1%	-7.5%	-13.2%	-20.1%	114k	35.0%
Scenario 3	326k	76.3%	-20.1%	-4.1%	-7.5%	-9.5%	0.0%	191k	58.8%
UT: Pine Valley
Scenario 1	153k	94.2%	-37.2%	-6.9%	-21.5%	-6.7%	0.0%	43k	27.8%
Scenario 2	153k	94.2%	-37.2%	-1.4%	-10.5%	-9.2%	0.0%	64k	41.7%
Scenario 3	153k	94.2%	-37.2%	-1.4%	-10.5%	-6.5%	0.0%	68k	44.4%
UT: Wasatch
Scenario 1	381k	89.5%	-43.0%	-12.1%	-14.8%	-4.9%	-0.3%	95k	25.0%
Scenario 2	381k	89.5%	-43.0%	-3.7%	-7.9%	-6.7%	-0.5%	146k	38.3%
Scenario 3	381k	89.5%	-43.0%	-3.7%	-7.9%	-4.8%	0.0%	155k	40.7%
WA: Central Washington Initiative
Scenario 1	893k	69.8%	-25.0%	-29.5%	-11.5%	-8.1%	-12.7%	117k	13.1%
Scenario 2	893k	69.8%	-25.0%	-9.5%	-7.0%	-13.4%	-22.2%	204k	22.9%
Scenario 3	893k	69.8%	-25.0%	-9.5%	-7.0%	-9.7%	0.0%	436k	48.8%
WA: Colville Northeast Washington Vision
Scenario 1	656k	92.8%	-20.1%	-18.5%	-16.0%	-10.1%	-1.2%	223k	34.1%
Scenario 2	656k	92.8%	-20.1%	-3.8%	-7.5%	-14.4%	-2.0%	342k	52.2%
Scenario 3	656k	92.8%	-20.1%	-3.8%	-7.5%	-10.4%	0.0%	382k	58.3%

4.4.2 Change in Available by Scenario

constrained_by_scnro_ls %>% 
  dplyr::mutate(scenario_lab = scenario_lab %>% forcats::fct_rev()) %>% 
ggplot(
    mapping = aes(x=scenario_lab,y=pct_rmn5_administrative,group=area_name)
    ) +
    geom_line(mapping=aes(color = area_name), linewidth = 1.5) +
    geom_label(
      mapping=aes(label = scales::percent(pct_rmn5_administrative, scale = 100, accuracy = 1))
      , color = "black"
      , size = 3
      , label.padding = unit(0.15, "lines")
    ) +
    facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , scales = "free_x"
    ) +
    scale_color_viridis_d(option = "turbo", alpha = 0.8) +
    scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)), labels = scales::percent_format(), breaks = scales::extended_breaks(6)) +
    # scale_x_discrete(labels = scales::label_wrap(20)) +
    labs(
      x = "Most Restrictive \U2192 Least Restrictive"
      , y = "Percent Treatable Area Remaining"
    ) +
    theme_light() +
    theme(
      legend.position = "none"
      , axis.text.y = element_text(size=7)
      , axis.text.x = element_text(size=7.5)
      , strip.text = element_text(color = "black", face = "bold", size = 10)
    )

4.4.3 Reduction by Constraint

# reshape
constrained_by_scnro_ls_long <- constrained_by_scnro_ls %>%
  dplyr::select(state, name, area_name, tidyselect::starts_with("scenario_"), tidyselect::starts_with("pct_rdctn")) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("pct_rdctn")
    , names_to = "constraint"
    , values_to = "pct_rdctn"
    , names_prefix = "pct_rdctn"
    , values_drop_na = F
  ) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
          , "total"   
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
          , "Total"
        )
      ) %>% forcats::fct_rev()
  ) %>% 
  dplyr::left_join(
    constrained_by_scnro_ls %>% dplyr::select(state,name,scenario_id,pct_rdctn_total)
    , by = join_by(state,name,scenario_id)
  )

# color pallete removing dark blue
n_temp <- ((constrained_by_scnro_ls_long$constraint %>% unique() %>% length())-1)*2
plasma_custom <- viridis::plasma(n = n_temp)[seq(2,n_temp,2)]
# scales::show_col(plasma_custom)
# plasma_custom
# plot
ggplot(
  data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint!="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
  , mapping = aes(y = area_name, x = pct_rdctn)
) +
  geom_col(
    mapping = aes(fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_rdctn<0.11*-1,pct_rdctn,NA), accuracy = 1)
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint=="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
    , mapping = aes(
      y = area_name, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(cols = vars(scenario_lab)) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.04, .3)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x =  element_blank() #
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_rdctn_ls.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

Adding in bar for remaining area:

constrained_by_scnro_ls_long %>% 
  dplyr::mutate(
    constraint = forcats::fct_recode(constraint,"Available"="Total")
    , pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(constraint)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = scenario_lab
    , x = pct_rdctn_hack
    , fill = constraint
    , group = constraint
  )
) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 1.8
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(constraint=="Available",pct_rdctn_hack,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = -0.1
  ) +
  annotate(
    geom = "text"
    , x = -0.78
    , y = 0
    , label = "constrained"
    , color = cols_cnstrnt_avail[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.88
    , y = 0
    , label = "available"
    , color = cols_cnstrnt_avail[2]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
  scale_x_continuous(
    limits = c(-1.00,1.02)
    , labels = scales::percent_format()
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = ""
    , x = "Percent Available and Constrained Relative to Total Forest and Shrub Area" 
    # , x = "Percent of Forest & Shrub Area in Overall Landscape" 
    , subtitle = "" 
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(fill = guide_legend(nrow = 1))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_avail.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

constrained_by_scnro_ls_long %>% 
  dplyr::mutate(
    constraint = forcats::fct_recode(constraint,"Available"="Total")
    , pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(constraint)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(scenario_lab, area_name)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = area_name
    , x = pct_rdctn_hack
    , fill = constraint
    , group = constraint
  )
) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "black", linetype = "solid", linewidth = 1.1) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(constraint=="Available",pct_rdctn_hack,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = -0.1
  ) +
  annotate(
    geom = "text"
    , x = -0.72
    , y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
    , label = "constrained"
    , color = cols_cnstrnt_avail[1]
    , fontface = "bold"
    , size = 3
  ) +
  annotate(
    geom = "text"
    , x = 0.82
    , y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
    , label = "available"
    , color = cols_cnstrnt_avail[2]
    , fontface = "bold"
    , size = 3
  ) +
  scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
  scale_x_continuous(
    limits = c(-1.00,1.02)
    , labels = scales::percent_format()
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.02,0.05))) +
  facet_grid(cols = vars(scenario_lab)) +
  labs(
    fill = ""
    , x = "Percent Available and Constrained Relative to Total Forest and Shrub Area" 
    # , x = "Percent of Forest & Shrub Area in Overall Landscape" 
    , subtitle = "" 
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(fill = guide_legend(nrow = 1))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_avail_wide.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 6
    , units = "in"
    , dpi = "print"
  )

4.5 Fireshed-level analysis

Based on model simulations of how much area generally needs to be treated to influence wildfire behavior, we binned the firesheds into three classes of mechanical constraint:

• high (81–100% [i.e., only 0–19% is available for mechanical treatment]): fuels treatment would principally need to rely on fire
• medium (60–80% [i.e., 20–40% is available]): could use a combination of fire and mechanical thinning
• low (60% [i.e., >40% is available]): could effectively influence wildfire behavior with mechanical treatment alone

4.5.1 Distribution of Fireshed PA Constraint

constrained_by_scnro_ls_pa %>% 
  dplyr::count(area_name,scenario_lab,cnstrnt_class) %>% 
  dplyr::group_by(area_name,scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , high_pct=max(ifelse(cnstrnt_class=="high constraint",pct,0)) 
    , tot = sum(n)
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.12,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  # geom_text(
  #     data = constrained_byftr_huc12_wide %>%
  #       dplyr::group_by(area_name) %>%
  #       dplyr::summarise(n=n(), high_n=sum(ifelse(cnstrnt_class=="high constraint",1,0)))
  #     , mapping = aes(x = -0.05,y=reorder(area_name, -(high_n/n)),label=paste0("n=",scales::comma(n)))
  #     , size = 3
  #     , color = "black"
  #     , hjust = 0.7
  #   ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(scenario_lab)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_fpa_cnstrnt_dist.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

4.5.2 Distribution of Fireshed Area

# plot
ggplot(data = constrained_by_scnro_ls_pa %>% dplyr::filter(scenario_id==1)
         , mapping = aes(x = feature_area_ha, y = reorder(area_name,desc(area_name) ))
    ) +
    geom_vline(
      xintercept = 
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::pull(feature_area_ha) %>% 
          median()
      , linetype="dashed", color="gray66"
    ) +
    geom_violin(aes(fill = state)) + 
    geom_boxplot(width = 0.15) +
    # geom_point(size = 0.7, color = cols_obj_r[2], alpha = 0.2) +
    geom_text(
      data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>% 
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(n=n(), max_area=max(feature_area_ha))
      , mapping = aes(
          x = (
            constrained_by_scnro_ls_pa %>% 
            dplyr::filter(scenario_id==1) %>% 
            dplyr::pull(feature_area_ha) %>% 
            min()
          )*.9
          ,y=reorder(area_name,desc(area_name) )
          ,label=paste0("n=",scales::comma(n))
        )
      , size = 3
      , color = "black"
      , hjust = 0.7
    ) +
    scale_fill_viridis_d(option = "viridis", alpha = 0.8) +
    scale_x_continuous(labels = scales::comma, breaks = scales::extended_breaks(n=7)) +
    labs(
      fill = ""
      , y = ""
      , x = "Fireshed Area (ha)"
    ) +
    theme_light() +
    theme(
      legend.position = "none" # c(0.9, 0.9)
      , axis.title = element_text(size=9)
      , axis.text.x = element_text(size=7)
    )

4.5.3 Map of Firesheds by Level of Constraint

# polish for mapping
map_firesheds <- fireshed_proj_area %>% 
  dplyr::select(pa_id, geometry) %>% 
  dplyr::mutate(pa_id=as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa
    , by = dplyr::join_by("pa_id")
    , multiple = "all"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(as.numeric(.x), accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::ends_with("_ha")
      , ~ scales::comma(as.numeric(.x), accuracy = 1)
    )
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
  dplyr::mutate(
    Priority.Area.Name = area_name
    , Level.of.Constraint = cnstrnt_class
    , Pct.Treatable.Class = rmn_cnstrnt_class
    , Fireshed.ProjArea.ID = pa_id
    , Fireshed.ProjArea.Area.Hectares = feature_area_ha
    , Fireshed.WCS.Status = ifelse(is.na(fireshed_crisis_strategy),"Not High Risk",fireshed_crisis_strategy)
    , ForestShrub.Area.Hectares = covertype_area_ha
    , Pct.Rdctn.Protected = pct_rdctn1_protected
    , Pct.Rdctn.Slope = pct_rdctn2_slope
    , Pct.Rdctn.Roads = pct_rdctn3_roads
    , Pct.Rdctn.Riparian = pct_rdctn4_riparian
    , Pct.Rdctn.Administrative = pct_rdctn5_administrative
    , Treatable.Area.Hectares = rmn5_administrative_area_ha
    , Pct.Treatable.Area = pct_rmn5_administrative
  )
# filter for mapview
mapview_fs_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 2
  )
# ## filter for high priority only
# hi_map_firesheds <- map_firesheds %>% 
#   dplyr::filter(!is.na(fireshed_crisis_strategy))
# filter for mapview
mapview_fs_sc1_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 1
  )
mapview_fs_sc3_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 3
  )
mapview_fs_sc3_hi_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 3
    & fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
  )

# mapview
mapview::mapviewOptions(homebutton = FALSE, basemaps = c("OpenStreetMap","Esri.WorldImagery"))
mapview::mapview(wf_landscapes
        , color = "black"
        , lwd = 1
        , alpha.regions = 0
        , label = FALSE
        , legend = FALSE
        , popup = FALSE
        , layer.name = "WCS Landscapes"
        , hide = TRUE
  ) +
mapview::mapview(
  mapview_fs_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 2 Constraints"
  , hide = FALSE
  , popup = leafpop::popupTable(
      mapview_fs_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc1_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 1 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc1_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc3_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) +
mapview::mapview(
  mapview_fs_sc3_hi_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints (High Only)"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_hi_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) 

# **Scenario 2 used for map above**

4.5.4 Map of Fireshed Project Areas by Scenario

plt_fn_temp <- function(scnum=1) {
  plt_basemap_simple +
    geom_sf(data = map_firesheds %>% dplyr::filter(scenario_id==scnum)
            , aes(fill = cnstrnt_class), lwd = 0.05) +
    
  geom_sf(
    data = wf_landscapes
    # , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    , color="black"
    , lwd=0.2
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  # scale_color_manual(values = c("black", cols_obj_r[2])) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  labs(
    fill = ""
    , tag = paste0(
      map_firesheds %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    # , plot.tag.position = c(0.8,0.94) # "top"
  ) + 
  guides(color = "none")
}
# plot
# plt_fn_temp()
plt_fpas_sc_map <- 
  map_firesheds %>% dplyr::pull(scenario_id) %>% unique() %>% 
    purrr::map(plt_fn_temp)
# print and export maps
  1:length(plt_fpas_sc_map) %>% 
    purrr::map(function(x){
    # save export images for publication
    print(plt_fpas_sc_map[x])
        ggplot2::ggsave(
            filename = paste0("../data/","plt_map_pa_cnstrnt_",x,".jpeg")
            # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
            , plot = ggplot2::last_plot()
            , width = 8
            , height = 8.5
            , units = "in"
            , dpi = "print"
            , bg = "white"
          )
      }
    )

4.5.5 Map Fireshed Project Areas by Scenario and Landscape

Define function to map each landscape separately with project areas highlighted by level of constraint. This function is used below to create figures with largest patch highlighted.

#########################################################################
#########################################################################
# Make each plot individually by landscape as solution to small multiples
# this block defines function
#########################################################################
##################hack to align plots for ggmap
ggmap_bbox_fn <- function(map, my_crs=3857) {
    if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
    # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
    # and set the names to what sf::st_bbox expects:
    map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
    # Convert the bbox to an sf polygon, transform it to 3857, 
    # and convert back to a bbox (convoluted, but it works)
    bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), my_crs))
    # Overwrite the bbox of the ggmap object with the transformed coordinates 
    attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
    attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
    attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
    attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
    map
}
plt_crs <- 3857
#########################################################################
#########################################################################
#########################################################################
area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
# buffer for smaller areas
if_smaller_ha_temp <- wf_landscapes %>% 
  dplyr::pull(hectares) %>% 
  quantile(0.45) %>% 
  purrr::pluck(1)
# create color palette
# function to plot basemap
landscape_basemap_fn <- function(row_n=1) {
  # should zoom in?
  zoom_level <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "WA: Colville Northeast Washington Vision"
      , "CA: Southern California Fireshed Risk Reduction Strategy"
      ) ~ 7
    , area_nm_list[row_n] %in% c(
      "UT: Wasatch"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "NV: Sierra and Elko Fronts"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 8
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "MT: Kootenai Complex"
      , "UT: Pine Valley"
    ) ~ 9
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 10
    , TRUE ~ 8
  )
  # should buffer extend?
  buffer_box <- dplyr::case_when(
    area_nm_list[row_n] %in% c("WA: Colville Northeast Washington Vision") ~ 58600
    , area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "MT: Kootenai Complex"
    ) ~ 40000
    , area_nm_list[row_n] %in% c(
      "NV: Sierra and Elko Fronts"
    ) ~ 21000
    , area_nm_list[row_n] %in% c(
      "NM: Enchanted Circle"
      , "UT: Pine Valley"
    ) ~ 25000
    , area_nm_list[row_n] %in% c(
      "AZ: San Carlos Apache Tribal Forest Protection"
      , "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 35500
    , area_nm_list[row_n] %in% c("UT: Wasatch") ~ 5000
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 15500
    , TRUE ~ 6000
  )
  
  # bounding box
  bb_temp <- 
    map_firesheds %>% 
    dplyr::filter(
      area_name == area_nm_list[row_n]
    ) %>% 
    sf::st_union() %>% 
    sf::st_transform(crs=5070) %>% 
    sf::st_buffer(as.numeric(buffer_box)) %>% 
    sf::st_transform(crs=4326) %>% # same as get_map return
    sf::st_bbox()
  # set bbox for get call
  bbox_temp <- c(
    bottom = bb_temp[[2]]
    , top = bb_temp[[4]]
    , right = bb_temp[[3]]
    , left = bb_temp[[1]]
  )
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # stamen maps is being depreciated Oct 31, 2023
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # get map
  hey_ggmap <- ggmap::get_stadiamap(
    bbox = bbox_temp
    , zoom = zoom_level
    , maptype = "stamen_terrain" #"toner-hybrid" # "terrain"
    , crop = T
  )

  # ggmap(hey_ggmap)
  # apply align function
    hey_ggmap_aligned <- ggmap_bbox_fn(hey_ggmap, plt_crs) # Use the function
  # plot
  plt <- ggmap(hey_ggmap_aligned) + 
    coord_sf(
      expand = FALSE
    ) +
    labs(
      title = area_nm_list[row_n]
    ) +
    theme_light() +
    theme(
      legend.position = "bottom"
      , legend.direction  = "horizontal"
      , legend.title = element_text(size=7)
      , legend.margin = margin(c(-10,0,0,0))
      , plot.title = element_text(size = 11)
      , strip.text = element_text(color = "black", face = "bold")
      , axis.title = element_blank()
      , axis.text = element_blank()
      , axis.ticks = element_blank()
      , panel.grid = element_blank()
      , plot.margin = margin(0, 0, 0, 0, "cm")
    )
  # return
  return(plt)
}
# landscape_basemap_fn(row_n = 13)
# function to plot scenario constraints + basemap
scenario_cnstrnt_map_fn <- function(row_n = 1) {
  # location ggspatial::annotation_scale 
  my_scale_loc <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "WA: Colville Northeast Washington Vision"
      , "CO: Colorado Front Range"
      , "AZ: 4FRI"
      , "ID: Southwest Idaho"
    ) ~ "bl"
    , TRUE ~ "br"
  )
  # plot
  plt <- landscape_basemap_fn(row_n) + 
      geom_sf(
        data = map_firesheds %>%
          dplyr::filter(
            area_name == area_nm_list[row_n]
          ) %>%
          sf::st_transform(crs=plt_crs)
        , aes(fill = cnstrnt_class), lwd = 0.05
        , inherit.aes = F
        , alpha = 0.6
      ) +
      geom_sf(
        data = wf_landscapes %>% 
          dplyr::filter(
              area_name == area_nm_list[row_n]
            ) %>% 
          sf::st_transform(crs=plt_crs)
        , fill = NA, color = "black", lwd = 0.3
        , inherit.aes = F
      ) +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      ggspatial::annotation_scale(
        location = my_scale_loc
        , style = "ticks"
        , pad_x = unit(0.1, "cm")
        , pad_y = unit(0.1, "cm")
      ) +
      labs(
        fill = "Level of\nConstraint"
      ) +
        guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(row_n = 15)

# same output with map of each landscape by scenario but without largest patch highlighted
# go through list of areas to make plot
plt_list <- 1:length(area_nm_list) %>%
  purrr::map(scenario_cnstrnt_map_fn)
# scenario_maps
plt_list[]
####
# area_nm_list
# scenario_cnstrnt_map_fn(18)

4.5.6 Patch Analysis

Notes from Collins et al. 2010:

The larger the individual fuel treatment the greater the potential for tree survival, because of reduced edge effect mortality (Ritchie et al. 2007). Larger individual treatments have a greater potential to reduce fire behavior and slow fire spread, which ultimately impacts adjacent untreated stands and should enhance suppression opportunities and increase firefighter safety (Collins et al. 2010). Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# create spatial data of interconnected patches by constraint class
# this is to represent the spatial aggregation of fireshed project areas of similar levels of  mechanical constraint
  # aggregation = Tendency of patch or land-cover types to be spatially adjacent or in close proximity
cnstrnt_class_patches <-
  fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::select(pa_id, area_name, scenario_id, scenario_lab, cnstrnt_class)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
  # union vectors/polygons by constraint class into one multipolygon
    # this dissolves inner boundaries by constraint class
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    geometry = sf::st_union(geometry, by_feature = F)
  ) %>% 
  # separate a multipolygon geometry into several polygons objects after performing a st_union()
  sf::st_cast(to = "MULTIPOLYGON") %>%
  sf::st_cast(to = "POLYGON") %>% 
  # calculate patch area 
  dplyr::ungroup() %>% 
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::mutate(
    patch_area_ha = as.numeric(sf::st_area(geometry))/10000
    , patch_area_rank = dplyr::dense_rank(dplyr::desc(patch_area_ha))
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(scenario_id, area_name, cnstrnt_class, dplyr::desc(patch_area_ha)) %>% 
  # calculate proportion of total area in each patch
  dplyr::inner_join(
    # total area of fireshed project areas included for analysis
      # after union to remove internal boundaries
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::select(pa_id, area_name)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
      dplyr::group_by(area_name) %>% 
      dplyr::summarize(
        geometry = sf::st_union(geometry, by_feature = F)
      ) %>%
      dplyr::mutate(
        total_area_ha = as.numeric(sf::st_area(geometry))/10000
      ) %>% 
      sf::st_drop_geometry()
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    pct_patch_landscape = patch_area_ha / total_area_ha
  )

# function to highlight largest patch with fireshed project areas under
largest_patch_map_fn <- function(nrow = 1){
    plt <- scenario_cnstrnt_map_fn(nrow) +
      geom_sf(
          data = cnstrnt_class_patches %>%
            dplyr::filter(
              area_name == area_nm_list[nrow]
              & patch_area_rank == 1
              & cnstrnt_class != "med. constraint"
            ) %>%
            sf::st_transform(crs=plt_crs)
          , aes(color = cnstrnt_class, fill = cnstrnt_class)
          , lwd = 0.9
          , inherit.aes = F
          , alpha = 0
          # , show.legend = F
        ) +
      scale_color_manual(values = cols_rdylbu_dk, drop = F) +
      labs(
        color = "Level of\nConstraint"
        , subtitle = "largest patch of interconnected fireshed project areas by level of mechanical constraint" # largest patch of spatially adjacent...
      ) +
      theme(
        # plot.subtitle = element_text(size = 8)
        plot.subtitle = element_blank()
      ) +
      guides(
        color = "none"
        , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(4)
# largest_patch_map_fn(5)

# plot each landscape
plt_ptch_list <- 1:length(area_nm_list) %>%
  purrr::map(largest_patch_map_fn)
# scenario_maps
plt_ptch_list[]

# save export images for publication
1:length(plt_ptch_list) %>% 
  purrr::map(function(x){
      # get height vs width
      x_temp = abs(plt_ptch_list[[x]]$data$lon[1]-plt_ptch_list[[x]]$data$lon[4])
      y_temp = abs(plt_ptch_list[[x]]$data$lat[1]-plt_ptch_list[[x]]$data$lat[4])
      # save as long or wide
      ggplot2::ggsave(
          filename = paste0("../data/patch_maps/","plt_map_patches_fpa_",x,".jpeg")
          # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
          , plot = plt_ptch_list[[x]] # ggplot2::last_plot()
          , width = ifelse(y_temp/x_temp > 1.5, 8.5, 11) # 11
          , height = ifelse(y_temp/x_temp > 1.5, 11, 8.5) # 8.5
          , units = "in"
          , dpi = "print"
          , bg = "white"
        )
    }
  )

4.5.6.1 Example of how patches were done

name_temp <- "WA: Central Washington Initiative"
# plot project areas and landscape
plt_patchex_stp1 <- ggplot() + 
  geom_sf(
    data = 
      fireshed_proj_area %>% 
        dplyr::mutate(pa_id=as.character(pa_id)) %>% 
        dplyr::select(pa_id) %>% 
        dplyr::inner_join(
          constrained_by_scnro_ls_pa %>% 
            dplyr::filter(
              area_name == name_temp
              & scenario_id == 2
            )
          , by = dplyr::join_by(pa_id)
        )
    , aes(fill = cnstrnt_class), lwd = 0.05
    , inherit.aes = F
    , alpha = 0.6
  ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp1
# plot patches and landscape
plt_patchex_stp2 <-
  ggplot() + 
  geom_sf(
      data = cnstrnt_class_patches %>%
        dplyr::filter(
          area_name == name_temp
          & scenario_id == 2
          # & patch_area_rank == 1
          # & cnstrnt_class != "med. constraint"
        )
      , aes(fill = cnstrnt_class) #, color = cnstrnt_class)
      # , lwd = 1.4
      , alpha = 0.6
    ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_dk, drop = F) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp2
# selected largest patch
plt_patchex_stp3 <-
  ggplot() + 
  geom_sf(
      data = cnstrnt_class_patches %>%
        dplyr::filter(
          area_name == name_temp
          & scenario_id == 2
          & patch_area_rank == 1
          & cnstrnt_class != "med. constraint"
        )
      , aes(fill = cnstrnt_class, color = cnstrnt_class)
      , lwd = 1.4
      , alpha = 0.6
    ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_dk, drop = F) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    color = "none"
    , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp3
### combine
cowplot::plot_grid(
  plotlist = list(plt_patchex_stp1, plt_patchex_stp2, plt_patchex_stp3)
  , cols = 3
  , rows = 1
  , align = "h"
)

ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_example_methods.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 7
    , units = "in"
    , dpi = "print"
  )

Comparison of largest patch of lightly constrained fireshed project areas versus largest patch of highly constrained fireshed project areas.

Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# aggregate to constraint class level
agg_cnstrnt_class_patches <- cnstrnt_class_patches %>% 
  sf::st_drop_geometry() %>% 
  dplyr::group_by(area_name, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    max_patch_area_ha = max(patch_area_ha, na.rm = T)
    , mean_patch_area_ha = mean(patch_area_ha, na.rm = T)
    , max_pct_patch_landscape = max(pct_patch_landscape, na.rm = T)
    , sum_pct_patch_landscape = sum(pct_patch_landscape)
  ) %>% 
  dplyr::ungroup() %>%
  dplyr::mutate(
    max_pct_patch_landscape_hack = dplyr::case_when(
      cnstrnt_class=="high constraint" ~ -max_pct_patch_landscape
      , TRUE ~ max_pct_patch_landscape
    )
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  )
# combine with full data to get 0's
agg_cnstrnt_class_patches <- 
  # full data
  dplyr::cross_join(
    constrained_by_scnro_ls_pa %>% dplyr::select(area_name) %>% dplyr::distinct()
    , constrained_by_scnro_ls_pa %>% dplyr::select(scenario_lab) %>% dplyr::distinct()
  ) %>% 
  dplyr::cross_join(
    constrained_by_scnro_ls_pa %>% dplyr::select(cnstrnt_class) %>% dplyr::distinct()
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  ) %>% 
  # end full data
  dplyr::left_join(
    agg_cnstrnt_class_patches
    , by = dplyr::join_by(area_name, scenario_lab, cnstrnt_class)
  ) %>% 
  # fill with 0's
  dplyr::mutate(
    dplyr::across(
      .cols = tidyselect::where(is.numeric) & -c(area_name, scenario_lab, cnstrnt_class, max_pct_patch_landscape_hack)
      , .fn = function(x){ifelse(is.na(x),0,x)}
    )
    , max_pct_patch_landscape_hack = ifelse(is.na(max_pct_patch_landscape_hack),-0.001,max_pct_patch_landscape_hack)
  )
# plot
ggplot(
  data = agg_cnstrnt_class_patches
  , mapping = aes(
    y = scenario_lab
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  # geom_line(
  #   mapping = aes(x = forcats::fct_rev(scenario_lab), y = max_pct_patch_landscape, color = forcats::fct_rev(cnstrnt_class), group = forcats::fct_rev(cnstrnt_class))
  #   , linewidth = 1.2
  # ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.73
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.92
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.22)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none" # "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_index.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

plot by scenario wide for presentation

# median for plot
median_temp <- agg_cnstrnt_class_patches %>%
        dplyr::mutate(
          dplyr::across(
            c(scenario_lab, area_name)
            , ~ forcats::fct_rev(.x)
          )
        ) %>% 
        dplyr::group_by(scenario_lab,cnstrnt_class) %>% 
        dplyr::summarise(
          median_pct = median(max_pct_patch_landscape_hack,na.rm=T)
        ) %>% 
        dplyr::ungroup()
# plot
ggplot(
  data = agg_cnstrnt_class_patches %>%
    dplyr::mutate(
      dplyr::across(
        c(scenario_lab, area_name)
        , ~ forcats::fct_rev(.x)
      )
    )
  , mapping = aes(
    y = area_name
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  geom_vline(
    # median line
    data = median_temp
    , mapping = aes(xintercept=median_pct, group = cnstrnt_class, color = cnstrnt_class)
    , linetype = "dashed"
    , lwd = 0.5
  ) +
  geom_text(
    # median line
    data = median_temp
    , mapping = aes(
      x=median_pct
      , y=0
      , group = cnstrnt_class
      , color = cnstrnt_class
      , label = scales::percent(median_pct, accuracy = 1)
      , fontface = "bold"
    )
    , size = 2
    , hjust = 1
  ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.75
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2
  ) +
  annotate(
    geom = "text"
    , x = 0.95
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.22)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.07,0.02))) +
  facet_grid(
    cols = vars(scenario_lab)
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_index_wide.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 6
    , units = "in"
    , dpi = "print"
  )

Combine map with proportion of overall landscape in largest patch plot for publication

patch_map_bars_combo_fn <- function(area_num_temp = 20) {

  # plot large patch proportion of area horizontal to match map
  plt_patch_prop_temp <- 
    ggplot(
      data = agg_cnstrnt_class_patches %>%
        dplyr::filter(area_name==area_nm_list[area_num_temp]) %>% 
        dplyr::mutate(scenario_lab = forcats::fct_rev(scenario_lab))
      , mapping = aes(
        y = 1
        , x = max_pct_patch_landscape_hack
        , fill = cnstrnt_class
        , group = cnstrnt_class
      )
    ) +
      geom_col(
        width = 0.7, alpha = 0.9
      ) +
      geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack<0
            , scales::percent(max_pct_patch_landscape, accuracy = 1)
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 4
        , hjust = +1
      ) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack>=0
            , scales::percent(max_pct_patch_landscape, accuracy = 1)
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 4
        , hjust = 0
      ) +
      # label ha below pct
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack<0
            , paste0(
                scales::comma(max_patch_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
                , " ha"
            )
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 2.7
        , hjust = +1
        , vjust = 1.9
      ) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack>=0
            , paste0(
                scales::comma(max_patch_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
                , " ha"
            )
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 2.7
        , hjust = 0
        , vjust = 1.9
      ) +
      annotate(
        geom = "text"
        , x = -0.85
        , y = 0.0
        , label = "high constraint"
        , color = cols_rdylbu_lt[3]
        , fontface = "bold"
        , size = 2.5
      ) +
      annotate(
        geom = "text"
        , x = 0.85
        , y = 0.0
        , label = "low constraint"
        , color = cols_rdylbu_lt[1]
        , fontface = "bold"
        , size = 2.5
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      scale_x_continuous(
        limits = c(-1.22,1.22)
        , labels = scales::percent_format()
        , position = "bottom"
      ) +
      scale_y_discrete() +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      labs(
        fill = "Level of\nConstraint"
        , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
        , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
        , y = ""
      ) +
      theme_light() +
      theme(
        legend.position = "none" # "top"
        , legend.direction  = "horizontal"
        , legend.title = element_text(size=7)
        , axis.title.x = element_text(size=8, face = "bold")
        , axis.title.y = element_blank()
        , axis.text.x = element_blank()
        , axis.ticks.x = element_blank()
        , panel.grid.major = element_blank()
        , panel.grid.minor = element_blank()
        , strip.text = element_text(color = "black", face = "bold", size = 8)
        # , strip.text = element_blank()
      ) +
      guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # combine
  plt <- cowplot::plot_grid(
      plotlist = list(
        plt_ptch_list[[area_num_temp]] + 
          theme(
            legend.position = "top"
            , plot.title = element_blank()
            , plot.margin = margin(0.5,0.1,0,0.1, "cm")
          )
        # , plt_fshed_pa + 
        #   theme(
        #     plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
        #     , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        #   )
        , plt_patch_prop_temp +
          theme(
            plot.margin = margin(0,0.9,0.5,0.9, "cm")
          )
        #######################################
        # , cowplot::plot_grid(
        #   plotlist = list(
        #     NULL
        #     , plt_patch_prop_temp +
        #       theme(
        #         plot.margin = margin(0,0.2,0.5,0.2, "cm")
        #       )
        #     , NULL
        #   )
        #   , nrow = 1
        #   , rel_widths = c(0.1,1,0.1)
        # )
        ######################################
        # , plt_fshed_pa_incl + 
        #   theme(
        #     plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
        #     , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        #   )
      )
      , ncol = 1
      , axis = "lr"
      , align = "v"
      , rel_heights = c(1, 0.25)
      , labels = c("A","B")
      , label_size = 17
      , label_fontface = "bold"
      , label_colour = "gray11"
      # , label_y = 0.99
    ) +
    theme(plot.background = element_rect(fill = "white", color = NA))
  # return plot
  return(plt)
}
# these are different sizes so manually set size
c(20) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 7.7
        , height = 7
        , units = "in"
        , dpi = "print"
      )
  })
c(9) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 7
        , height = 7.3
        , units = "in"
        , dpi = "print"
      )
  })
c(6) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 9
        , height = 6.4
        , units = "in"
        , dpi = "print"
      )
  })
# c(20,9,6) %>% 
#   purrr::map(function(x){
#     # export
#     ggplot2::ggsave(
#         filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
#         , plot = patch_map_bars_combo_fn(x)
#         , width = 7.7
#         , height = 7
#         , units = "in"
#         , dpi = "print"
#       )
#   })

4.6 High-risk fireshed project areas

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the plan (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

4.6.1 Project Areas by WCS Status

# aggregate data by area, status
wcs_status_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(scenario_id==1) %>%
  dplyr::group_by(area_name,fireshed_crisis_strategy) %>% 
  dplyr::summarise(
    n = n()
    , feature_area_ha = sum(feature_area_ha)
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    tot_n = sum(n)
    , tot_area = sum(feature_area_ha)
    , pct_n = n/tot_n
    , pct_area = feature_area_ha/tot_area
  ) %>% 
  dplyr::ungroup()
# plot
ggplot(data = wcs_status_temp) +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_n), x = n, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = n, group = fireshed_crisis_strategy
      , label = scales::comma(ifelse(n>20,n,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,n), x = n
      , label = scales::comma(n, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.1)),labels = scales::comma_format()) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "# of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

4.6.2 Area of high-risk project areas

Only including fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed.

See this section for discussion of high-risk firesheds and criteria for including fireshed project areas.

# data filter
wcs_status_temp %>% 
  dplyr::filter(fireshed_crisis_strategy!="Not High Risk") %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(tot_area = sum(feature_area_ha)) %>% 
# plot
ggplot() +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_area), x = feature_area_ha, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = feature_area_ha, group = fireshed_crisis_strategy
      , label = ifelse(
          feature_area_ha<150000
          ,NA
          ,scales::comma(feature_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
        )
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(
          scenario_id==1
          & fireshed_crisis_strategy!="Not High Risk"
        ) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,feature_area_ha), x = feature_area_ha
      , label = scales::comma(feature_area_ha,suffix = "k ha", scale = 1e-3, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15)),labels = scales::comma_format(suffix = "k ha", scale = 1e-3, accuracy = 1)) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "Area (ha) of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

4.6.3 Reduction by Constraint

See this section for reduction by constraint at the landscape-level.

# aggregate by area_name, scenario for hi-risk pa's
constrained_by_scnro_ls_hionly_long <-
  constrained_by_scnro_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha
  ) %>%
  dplyr::group_by(
    scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
  ) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  )
  
# plot
ggplot(data = constrained_by_scnro_ls_hionly_long) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_hionly_long %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = -1*(1 - pct_rmn)
      , label = scales::percent(-1*(1 - pct_rmn), accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10) 
  )

The figure above shows the reduction by constraint for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

4.6.4 Reduction by Constraint high-risk vs. total

# combine data and plot
pct_rdctn_temp <- 
constrained_by_scnro_ls_hionly_long %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_long %>% 
      dplyr::ungroup() %>% 
      dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
      dplyr::filter(!is.na(constraint_lvl)) %>% 
      dplyr::mutate(dta_src = "Overall Landscape Area")
  ) %>% 
  dplyr::group_by(scenario_lab, area_name, dta_src) %>% 
  dplyr::mutate(
    pct_rdctn_total = sum(pct_rdctn)
    , constraint = factor(
        constraint_lvl
        , ordered = TRUE
        , levels = 1:5
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  )
# plot
ggplot(data = pct_rdctn_temp) +
  geom_col(
    mapping = aes(y = dta_src, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = dta_src, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.075*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = pct_rdctn_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = dta_src, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "y"
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.text.y.left = element_text(angle = 0)
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

4.7 Objective vs. Treatable

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

Using the area of fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed as the treatment area objective (20-40%).

# Area of pa's by level of constraint
constrained_by_scnro_ls_cnstrnt_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name,cnstrnt_class) %>% 
  dplyr::summarise(
    n = n()
    , dplyr::across(
      c(feature_area_ha,covertype_area_ha,rmn5_administrative_area_ha)
      , ~ sum(.x,na.rm = T)
    )
  ) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name) %>% 
  dplyr::mutate(
    # cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
    total_area_ha = sum(feature_area_ha)
    , total_covertype_area_ha = sum(covertype_area_ha)
    # "the need to treat 20 to 40 percent of the overall fireshed"
    , pct_treatable =  rmn5_administrative_area_ha/total_area_ha
    , objective_lo_treat_area_ha = total_area_ha*.2
    , objective_hi_treat_area_ha = total_area_ha*.4
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(state,name,scenario_id,cnstrnt_class)

4.7.1 Objective vs. Treatable Table

table_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::mutate(
    cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
  ) %>% 
  dplyr::rename(
    remaining_pct=pct_treatable
    , remaining_ha=rmn5_administrative_area_ha
  ) %>%
  dplyr::select(-c(feature_area_ha, covertype_area_ha)) %>% 
  tidyr::pivot_wider(
    names_from = cnstrnt_class
    , values_from = c(n, remaining_pct, remaining_ha)
    , names_glue = "{cnstrnt_class}_{.value}"
    , values_fill = 0
  ) %>% 
  dplyr::mutate(
    tot_remaining_ha = low_remaining_ha+med_remaining_ha
    ,tot_remaining_pct = low_remaining_pct+med_remaining_pct
  ) %>% 
  dplyr::relocate(tot_remaining_ha,.after = med_remaining_ha) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name,hectares) %>% 
      dplyr::rename(landscape_area_ha=hectares)
    , by = dplyr::join_by(area_name)
  )
# make table 
kableExtra::kable(
    table_temp %>% 
      dplyr::select(
        c(
          scenario_lab
          # , landscape_area_ha
          , total_area_ha
          , total_covertype_area_ha
          # light and mod area
          , low_remaining_ha, med_remaining_ha, tot_remaining_ha
          , low_remaining_pct, med_remaining_pct, tot_remaining_pct
        )
      ) %>% 
      dplyr::mutate(
        dplyr::across(
          tidyselect::ends_with("_ha")
          , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
        )
        , dplyr::across(
          tidyselect::ends_with("_pct")
          , ~ scales::percent(.x, accuracy = .1)
        )
      )
  , caption = "Wildfire Crisis Strategy High-Risk Fireshed Project Areas<br>mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)"
  , col.names = c(
    ""
    # , "Total\nLandscape"
    , "Total Area (ha)"
    , "Forest+Shrub Area (ha)"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
  
  )
  , escape = F
) %>%  
    add_header_above(c(" " = 1, "High-Risk Fireshed\nProject Areas" = 2, "Treatable High-Risk\nArea (ha)"=3, "Treatable High-Risk\nPercent (%)" = 3)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(table_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::column_spec(
      9 # 10
      , color = "white"
      , background = 
          table_temp %>% 
            dplyr::mutate(
              is_gt20 = dplyr::case_when(
                low_remaining_pct >= .2 ~ cols_obj_html[1]
                , tot_remaining_pct >= .2 ~ cols_obj_html[2]
                , TRUE ~ cols_obj_html[3]
              )
          ) %>% 
          dplyr::pull(is_gt20)
          
    ) %>% 
    kableExtra::scroll_box(width = "740px")

Table 4.5: Table 4.6: Wildfire Crisis Strategy High-Risk Fireshed Project Areas
mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)

	High-Risk Fireshed Project Areas		Treatable High-Risk Area (ha)			Treatable High-Risk Percent (%)
	Total Area (ha)	Forest+Shrub Area (ha)	Lightly Constrained	Moderately Constrained	Total	Lightly Constrained	Moderately Constrained	Total
AZ: 4FRI
Scenario 1	1,579k	1,436k	304k	126k	430k	19.3%	8.0%	27.2%
Scenario 2	1,579k	1,436k	535k	69k	605k	33.9%	4.4%	38.3%
Scenario 3	1,579k	1,436k	913k	50k	963k	57.8%	3.2%	61.0%
AZ: Prescott
Scenario 1	300k	263k	34k	30k	63k	11.3%	9.9%	21.1%
Scenario 2	300k	263k	90k	25k	115k	29.9%	8.5%	38.4%
Scenario 3	300k	263k	134k	15k	150k	44.8%	5.0%	49.8%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1	250k	230k	4k	15k	19k	1.6%	6.1%	7.7%
Scenario 2	250k	230k	28k	50k	78k	11.3%	19.9%	31.2%
Scenario 3	250k	230k	69k	31k	100k	27.6%	12.2%	39.8%
CA: North Yuba
Scenario 1	206k	197k	21k	36k	57k	10.4%	17.3%	27.6%
Scenario 2	206k	197k	75k	13k	88k	36.5%	6.5%	43.0%
Scenario 3	206k	197k	99k	12k	111k	48.4%	5.7%	54.1%
CA: Plumas Community Protection
Scenario 1	157k	148k	12k	34k	46k	7.8%	21.7%	29.5%
Scenario 2	157k	148k	65k	8k	72k	41.2%	4.9%	46.1%
Scenario 3	157k	148k	91k	0k	91k	57.9%	0.0%	57.9%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1	1,053k	854k	13k	92k	105k	1.3%	8.7%	10.0%
Scenario 2	1,053k	854k	162k	87k	249k	15.3%	8.3%	23.7%
Scenario 3	1,053k	854k	241k	71k	311k	22.8%	6.7%	29.6%
CA: Stanislaus
Scenario 1	153k	143k	20k	30k	50k	13.1%	19.8%	32.8%
Scenario 2	153k	143k	59k	10k	70k	38.8%	6.7%	45.5%
Scenario 3	153k	143k	80k	4k	84k	52.4%	2.5%	54.9%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1	565k	488k	3k	32k	36k	0.6%	5.7%	6.3%
Scenario 2	565k	488k	49k	48k	97k	8.6%	8.5%	17.1%
Scenario 3	565k	488k	171k	48k	219k	30.3%	8.5%	38.8%
CO: Colorado Front Range
Scenario 1	1,138k	823k	85k	118k	203k	7.5%	10.4%	17.8%
Scenario 2	1,138k	823k	286k	65k	351k	25.1%	5.7%	30.8%
Scenario 3	1,138k	823k	374k	62k	436k	32.9%	5.4%	38.3%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1	520k	444k	77k	49k	126k	14.8%	9.5%	24.2%
Scenario 2	520k	444k	199k	13k	212k	38.3%	2.5%	40.8%
Scenario 3	520k	444k	205k	14k	218k	39.4%	2.6%	42.0%
ID: Southwest Idaho
Scenario 1	734k	632k	228k	41k	268k	31.0%	5.5%	36.5%
Scenario 2	734k	632k	408k	15k	423k	55.6%	2.1%	57.6%
Scenario 3	734k	632k	427k	8k	435k	58.1%	1.1%	59.2%
MT: Kootenai Complex
Scenario 1	322k	294k	98k	12k	111k	30.6%	3.9%	34.4%
Scenario 2	322k	294k	140k	9k	149k	43.5%	2.9%	46.4%
Scenario 3	322k	294k	190k	9k	199k	59.1%	2.7%	61.8%
NV: Sierra and Elko Fronts
Scenario 1	678k	453k	99k	55k	154k	14.6%	8.2%	22.8%
Scenario 2	678k	453k	229k	17k	247k	33.8%	2.6%	36.4%
Scenario 3	678k	453k	250k	20k	270k	36.9%	2.9%	39.9%
NM: Enchanted Circle
Scenario 1	318k	285k	33k	42k	75k	10.4%	13.1%	23.5%
Scenario 2	318k	285k	118k	16k	134k	37.1%	4.9%	42.0%
Scenario 3	318k	285k	133k	13k	146k	41.9%	4.0%	45.9%
OR: Central Oregon
Scenario 1	278k	232k	138k	10k	149k	49.7%	3.6%	53.4%
Scenario 2	278k	232k	169k	0k	169k	60.5%	0.0%	60.5%
Scenario 3	278k	232k	203k	0k	203k	73.1%	0.0%	73.1%
OR: Klamath River Basin
Scenario 1	684k	557k	13k	45k	59k	2.0%	6.6%	8.6%
Scenario 2	684k	557k	137k	36k	173k	20.0%	5.2%	25.3%
Scenario 3	684k	557k	220k	35k	254k	32.1%	5.1%	37.2%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1	339k	244k	15k	22k	37k	4.5%	6.5%	11.0%
Scenario 2	339k	244k	46k	21k	67k	13.6%	6.3%	19.9%
Scenario 3	339k	244k	131k	10k	141k	38.6%	2.8%	41.5%
UT: Pine Valley
Scenario 1	187k	173k	7k	32k	39k	3.9%	16.9%	20.8%
Scenario 2	187k	173k	36k	24k	60k	19.1%	13.1%	32.2%
Scenario 3	187k	173k	45k	19k	64k	24.0%	10.3%	34.4%
UT: Wasatch
Scenario 1	156k	126k	21k	11k	32k	13.6%	7.1%	20.7%
Scenario 2	156k	126k	45k	8k	52k	28.6%	5.0%	33.7%
Scenario 3	156k	126k	47k	8k	56k	30.4%	5.3%	35.7%
WA: Central Washington Initiative
Scenario 1	903k	654k	10k	37k	47k	1.1%	4.1%	5.2%
Scenario 2	903k	654k	73k	46k	119k	8.1%	5.1%	13.2%
Scenario 3	903k	654k	285k	43k	328k	31.6%	4.7%	36.3%
WA: Colville Northeast Washington Vision
Scenario 1	272k	244k	55k	28k	83k	20.4%	10.2%	30.6%
Scenario 2	272k	244k	118k	10k	128k	43.3%	3.6%	46.9%
Scenario 3	272k	244k	149k	6k	155k	54.7%	2.2%	56.9%

4.7.2 Reduction in treatable area flow

table_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::select(
    c(
      area_name
      , landscape_area_ha
      , total_area_ha
      , total_covertype_area_ha
      # , tot_remaining_ha
      # , low_remaining_ha
    )
  ) %>%
  tidyr::pivot_longer(
    !area_name
    , names_to = "name"
    , values_to = "hectares"
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(order = dplyr::row_number(),y_val=2) %>% 
  dplyr::ungroup() %>% 
  dplyr::cross_join(data.frame(scenario_id = 1:3)) %>% 
  ####
  dplyr::bind_rows(
    table_temp %>% 
      dplyr::select(
        c(
          area_name
          , scenario_id
          , tot_remaining_ha
          , low_remaining_ha
        )
      ) %>%
      tidyr::pivot_longer(
        !c(area_name,scenario_id)
        , names_to = "name"
        , values_to = "hectares"
      ) %>% 
      dplyr::group_by(area_name,scenario_id) %>% 
      dplyr::mutate(
        order = dplyr::row_number()+3
        , y_val = scenario_id
      ) %>% 
      dplyr::ungroup() %>% 
      dplyr::relocate(scenario_id,.after = dplyr::last_col())
  ) %>% 
  dplyr::arrange(area_name,order,y_val ) %>% 
  dplyr::mutate(
    order_lab = factor(
      order
      , levels = 1:5
      , labels = c(
        "Total\nLandscape"
        , "High-risk\nFireshed PAs"
        , "High-risk\nForest+Shrub"
        , "Lt.+Md. Constrained"
        , "Lightly Constrained"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name, hectares) %>% 
      dplyr::rename(tot_hectares=hectares)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::inner_join(
    table_temp %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(area_name, objective_lo_treat_area_ha)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    y_val = -y_val
    , sizing = ifelse(hectares>tot_hectares,1,hectares/tot_hectares)
    , coloring = dplyr::case_when(
        order %in% 4:5
          & hectares >= objective_lo_treat_area_ha
          ~ 1
        , order %in% 4:5
          & hectares < objective_lo_treat_area_ha
          ~ 2
        , T ~ 3
    )
  ) %>% 
  # dplyr::filter(stringr::str_starts(area_name,"CA")) %>% 
  # View()
ggplot() +
  geom_line(
    mapping = aes(x = order_lab, y = y_val, group = scenario_id)
    , color = "gray"
  ) +
  geom_point(
    mapping = aes(x = order_lab, y = y_val, size = sizing*100, color = factor(coloring))
    , shape = 15
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = scales::comma(hectares,suffix = "k", scale = 1e-3, accuracy = 1)
    )
    , size = 2.5
    , color = "black"
    , vjust = 2.5
    , angle = 90
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = ifelse(
        order == 4, paste0("Scenario", y_val*-1), NA
      )
      , fontface = "italic"
    )
    , size = 2.5
    , color = "black"
    , hjust = 0.1
    , vjust = -1.4
  ) +
  facet_wrap(facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_x"
  ) +
  scale_color_manual(values = c(cols_obj_html[1],cols_obj_html[3],cols_obj_html[2])) +
  scale_y_continuous(expand = expansion(c(0.3,0.3))) +
  scale_x_discrete(labels = scales::label_wrap(10)) +
  labs(
    x = ""
    , y = ""
    , subtitle = "<span><span style='color:#008b8b;'><b><i>≥20% treatable</i></b></span> | <span style='color:#cd3700;'><b><i><20% treatable</i></b></span></span>"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.y = element_blank()
    , axis.ticks.y = element_blank()
    , axis.text = element_text(size = 7)
    , panel.grid = element_blank()
    , plot.subtitle = ggtext::element_markdown(size = 10)
    , strip.text = element_text(color = "black", face = "bold", size = 10)
  )

4.7.3 Objective vs. Treatable Plot

obj_trtbl_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::filter(cnstrnt_class!="high constraint") %>% 
  dplyr::left_join(
    # aggregate by area_name, scenario_lab
    constrained_by_scnro_ls_cnstrnt_temp %>% 
        dplyr::filter(cnstrnt_class!="high constraint") %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::mutate(low=ifelse(cnstrnt_class=="low constraint",pct_treatable,0)) %>% 
        dplyr::summarise(
          pct_treatable_lowmed = sum(pct_treatable)
          , pct_treatable_low = sum(low)
        ) %>% 
        dplyr::mutate(
          treatment_objective = dplyr::case_when(
              pct_treatable_low >= .2 ~ 1
              , pct_treatable_lowmed >= .2 ~ 2
              , TRUE ~ 3
            ) %>% 
            factor(
              levels = 1:3
              , labels = c(
                  ">20% in\nlow constraint"
                  ,">20% in\nlow+medium constraint"
                  , "<20% in\nlow+medium constraint"
              )
              , ordered = T
            )
        ) %>% 
        dplyr::group_by(scenario_id) %>% 
        dplyr::arrange(scenario_id, treatment_objective
                       , desc(pct_treatable_low), desc(pct_treatable_lowmed)
        ) %>% 
        dplyr::mutate(sorting_rank = dplyr::row_number()) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(area_name, scenario_id)
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(
    dplyr::across(
      c(cnstrnt_class, scenario_lab)
      , ~ forcats::fct_rev(.x)
    )
  ) 
# plot
obj_plt_fn <- function(scnum){
plt <- ggplot(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum)
    , mapping = aes(
        x = pct_treatable
        , y = reorder(area_name, desc(sorting_rank))
      )
  ) +
  geom_rect(mapping = aes(xmin = 0.2, xmax=0.4, ymin = -Inf, ymax = Inf)
    , fill = "gray77"
    , alpha = 0.6
  ) +
  annotate("text"
     , x = 0.46
     , y = obj_trtbl_temp %>% 
          dplyr::filter(scenario_id == scnum & sorting_rank==max(obj_trtbl_temp$sorting_rank)) %>% 
          dplyr::pull(area_name) %>% 
          unique()
     , label = expression(bold("objective")~bold("20-40%"))
     , size = 3, color = "gray65", parse = T
    ) +
  geom_col(
    mapping = aes(fill = cnstrnt_class)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      group = cnstrnt_class
      , label = scales::percent(
          ifelse(pct_treatable<.06,NA,pct_treatable)
          , accuracy = 1
        )
      # , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.2
  ) +
  # total
  geom_text(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
    , mapping = aes(
      x = pct_treatable_lowmed
      , y = reorder(area_name, desc(sorting_rank))
      , label = scales::percent(pct_treatable_lowmed, accuracy = 1) 
      , color = treatment_objective
      , fontface = "bold"
    )
    , size = 3.5
    , hjust = -0.1
    , show.legend = F
  ) +
  # facet_grid(
  #   rows = vars(scenario_lab)
  #   , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
  #   , scales = "free_y"
  #   # , switch = "y"
  # ) +
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=3,name="RdYlBu")[2:3]) +
  scale_color_manual(values = c(
    ">20% in\nlow constraint"=cols_obj_r[1]
    ,">20% in\nlow+medium constraint"=cols_obj_r[2]
    , "<20% in\nlow+medium constraint"=cols_obj_r[3])
  ) +
  # scale_x_continuous(expand = expansion(mult = c(0, .15)),labels = scales::percent_format()) +
  scale_x_continuous(limits = c(0, .78),labels = scales::percent_format(),position = "top") +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Total High-Risk Area"
    , subtitle = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=9, face = "bold")
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , axis.text = element_text(size = 8)
    , axis.text.y = element_text(color = "black",size=9, face = "bold")
    , panel.grid.major.x = element_blank()
    , panel.grid.minor.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10)
    , strip.background = element_rect(fill = "gray88")
    # , strip.placement = "outside"
    , plot.margin = margin(0, 0, 0, 0, "cm")
    , plot.subtitle = element_text(vjust = -4,face="bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
return(plt)
}
cowplot::plot_grid(
  obj_plt_fn(1) +
    theme(
      legend.position = c(0.86,0.8)
      , legend.direction = "vertical"
    )
  , obj_plt_fn(2) +
    theme(legend.position = "none")
  , obj_plt_fn(3) +
    # scale_x_continuous(position = "top") +
    # labs(x = "Percent of Total High-Risk Area") +
    theme(legend.position = "none", plot.margin = margin(0, 0, 0.1, 0, "cm"))
  , ncol = 1
  , align = 'vh'
)

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_obj_sc.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

map landscapes shaded by potential to meet 20-40%

# map landscapes shaded by potential to meet 20-40%
plt_fn_temp <- function(scnum=1) {
  # filter and join data
  sf_dta_temp <- wf_landscapes %>% 
  dplyr::select(area_name) %>% 
    # join with spatial data
    dplyr::inner_join(
      obj_trtbl_temp %>% 
        dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by("area_name")
    ) 
  # plot
  plt <- plt_basemap_simple +
  geom_sf(
    data = sf_dta_temp
    , mapping = aes(fill = treatment_objective)
    , color="black"
    , lwd=0.1
    , alpha = 0.9
  ) +
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c(
    ">20% in\nlow constraint"=cols_obj_r[1]
    ,">20% in\nlow+medium constraint"=cols_obj_r[2]
    , "<20% in\nlow+medium constraint"=cols_obj_r[3])
  ) +
  labs(
    fill = ""
    , tag = paste0(
      map_firesheds %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    , legend.spacing.y = unit(0.25, 'cm')
    , plot.margin = margin(0,0,0,0,"cm")
  ) +
  guides(fill = guide_legend(byrow = TRUE))
  # return
  return(plt)
}
# plot
# save export images for publication
obj_trtbl_temp$scenario_id %>%
  unique() %>%
  purrr::map(function(x){
      print(plt_fn_temp(x))
      ggplot2::ggsave(
          filename = paste0("../data/","plt_map_objctv_ls_",x,".jpeg")
          # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
          , plot = ggplot2::last_plot()
          , width = 8
          , height = 8.5
          , units = "in"
          , dpi = "print"
          , bg = "white"
        )
    }
  )

make table of landscapes shifting availability toward objective

# kable
tbl_temp <- obj_trtbl_temp %>% 
  dplyr::group_by(area_name, scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::group_by(scenario_id, treatment_objective) %>% 
  # dplyr::arrange(desc(pct_treatable_low)) %>% 
  dplyr::arrange(area_name) %>% 
  dplyr::mutate(rnk = dplyr::row_number()) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, treatment_objective,rnk) %>% 
  tidyr::pivot_wider(
    names_from = scenario_lab
    , values_from = area_name
  ) %>% 
  dplyr::arrange(treatment_objective, rnk)
# kable table
  options(knitr.kable.NA = "") 
  kableExtra::kable(
    tbl_temp %>% dplyr::select(tidyselect::starts_with("Scenario"))
    , caption = " "
    # , col.names = c("")
  ) %>% 
  kable_classic("striped", full_width=T) %>%
  pack_rows(index = table(forcats::fct_inorder(tbl_temp$treatment_objective)))

Table 4.7: Table 4.8:

Scenario 1	Scenario 2	Scenario 3
>20% in low constraint
ID: Southwest Idaho	AZ: 4FRI	AZ: 4FRI
MT: Kootenai Complex	AZ: Prescott	AZ: Prescott
OR: Central Oregon	CA: North Yuba	AZ: San Carlos Apache Tribal Forest Protection
WA: Colville Northeast Washington Vision	CA: Plumas Community Protection	CA: North Yuba
	CA: Stanislaus	CA: Plumas Community Protection
	CO: Colorado Front Range	CA: Southern California Fireshed Risk Reduction Strategy
	ID: Nez Perce-Clearwater-Lower Salmon	CA: Stanislaus
	ID: Southwest Idaho	CA: Trinity Forest Health and Fire-Resilient Rural Communities
	MT: Kootenai Complex	CO: Colorado Front Range
	NM: Enchanted Circle	ID: Nez Perce-Clearwater-Lower Salmon
	NV: Sierra and Elko Fronts	ID: Southwest Idaho
	OR: Central Oregon	MT: Kootenai Complex
	OR: Klamath River Basin	NM: Enchanted Circle
	UT: Wasatch	NV: Sierra and Elko Fronts
	WA: Colville Northeast Washington Vision	OR: Central Oregon
		OR: Klamath River Basin
		OR: Mount Hood Forest Health and Fire-Resilient Communities
		UT: Pine Valley
		UT: Wasatch
		WA: Central Washington Initiative
		WA: Colville Northeast Washington Vision
>20% in low+medium constraint
AZ: 4FRI	AZ: San Carlos Apache Tribal Forest Protection
AZ: Prescott	CA: Southern California Fireshed Risk Reduction Strategy
CA: North Yuba	UT: Pine Valley
CA: Plumas Community Protection
CA: Stanislaus
ID: Nez Perce-Clearwater-Lower Salmon
NM: Enchanted Circle
NV: Sierra and Elko Fronts
UT: Pine Valley
UT: Wasatch
<20% in low+medium constraint
AZ: San Carlos Apache Tribal Forest Protection	CA: Trinity Forest Health and Fire-Resilient Rural Communities
CA: Southern California Fireshed Risk Reduction Strategy	OR: Mount Hood Forest Health and Fire-Resilient Communities
CA: Trinity Forest Health and Fire-Resilient Rural Communities	WA: Central Washington Initiative
CO: Colorado Front Range
OR: Klamath River Basin
OR: Mount Hood Forest Health and Fire-Resilient Communities
WA: Central Washington Initiative

4.8 Totals

##################
# Data Prep
##################
# union landscape and pa level data
# aggregate by scenario for hi-risk pa's
pct_rdctn_pas_temp <- 
  constrained_by_scnro_ls_pa %>%
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# overall
pct_rdctn_ls_temp <- 
  constrained_by_scnro_ls %>%
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# combine
pct_rdctn_all_temp <- 
  pct_rdctn_pas_temp %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
    dplyr::bind_rows(
      pct_rdctn_ls_temp %>% 
        dplyr::mutate(dta_src = "Overall Landscape Area")
    ) %>% 
    dplyr::mutate(
      constraint_lab = factor(
          constraint_lvl
          , ordered = TRUE
          , 1:5
          , labels = c(
            "Protected"
            , "Slope\nSteepness"
            , "Road\nDistance"
            , "Riparian\nBuffer"
            , "Administrative"
          )
        ) %>% forcats::fct_rev()
    ) %>% 
    dplyr::relocate(dta_src)

4.8.1 Reduction Treatable Area Table

tot_table_temp <- pct_rdctn_all_temp %>% 
  dplyr::mutate(constraint_nm = paste0(constraint_lvl,constraint)) %>% 
  dplyr::select(-c(constraint_lvl,constraint,constraint_lab)) %>% 
  tidyr::pivot_wider(
    names_from = constraint_nm
    , values_from = c(rmn_area_ha, pct_rmn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    pct_covertype_area = covertype_area_ha/feature_area_ha
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::contains("area_ha")
      , ~ scales::comma(.x,suffix = " M", scale = 1e-6, accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(.x, accuracy = 0.1)
    )
  ) %>% 
  dplyr::arrange(
    desc(dta_src)
    , scenario_id
  )
# make table 
  kableExtra::kable(
      tot_table_temp %>% 
        dplyr::select(
          scenario_lab
          , covertype_area_ha
          , pct_covertype_area
          , pct_rdctn_1protected
          , pct_rdctn_2slope
          , pct_rdctn_3roads
          , pct_rdctn_4riparian
          , pct_rdctn_5administrative
          , rmn_area_ha_5administrative
          , pct_rmn_5administrative
        )
      , caption = paste0("<b>Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas</b><br>percent reduction of different types of constraints on mechanical treatment")
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Strategy Priority\nArea"=2, "Constraints on\nMechanical Treatment" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tot_table_temp$dta_src))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")

Table 4.9: Table 4.10: Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas
percent reduction of different types of constraints on mechanical treatment

	Strategy Priority Area		Constraints on Mechanical Treatment
	Forest+Shrub (ha)	|Forest+Shrub	|Protecte	|Slope Steepne	s |Road Dista	ce |Riparian Bu	fer |Administra	ive |Remaining	ha) |Remaining
Overall Landscape Area
Scenario 1	16.1 M	82.9%	-22.0%	-16.9%	-19.8%	-7.0%	-5.6%	4.6 M	28.7%
Scenario 2	16.1 M	82.9%	-22.0%	-5.2%	-12.0%	-10.5%	-8.4%	6.8 M	41.9%
Scenario 3	16.1 M	82.9%	-22.0%	-5.2%	-12.0%	-7.5%	0.0%	8.6 M	53.3%
High-Risk Project Areas
Scenario 1	8.8 M	82.5%	-17.4%	-21.0%	-18.2%	-8.5%	-6.3%	2.5 M	28.5%
Scenario 2	8.8 M	82.5%	-17.4%	-6.4%	-10.5%	-13.0%	-9.5%	3.8 M	43.1%
Scenario 3	8.8 M	82.5%	-17.4%	-6.4%	-10.5%	-9.3%	0.0%	5.0 M	56.4%

4.8.2 Reduction by Constraint

# plot
# plot
ggplot(data = pct_rdctn_all_temp) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint_lab)
    , color = NA, width = 0.6
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = pct_rdctn_all_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(dta_src %>% forcats::fct_rev())
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0, .12)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

4.8.3 Distribution of Fireshed PA Constraint

cnstrnt_dta_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(scenario_desc, scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , dta_src = "High-Risk Firesheds"
  ) %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_pa %>% 
      dplyr::group_by(pa_id,scenario_id) %>% 
      dplyr::filter(dplyr::row_number()==1) %>% 
      dplyr::ungroup() %>% 
      dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
      dplyr::group_by(scenario_desc, scenario_lab) %>% 
      dplyr::mutate(
        pct = n/sum(n)
        , dta_src = "Overall Landscape Area"
      )
  ) %>% 
  dplyr::mutate(
    dta_src = dta_src %>% forcats::fct_rev()
  )
# n for notation
nlab_temp <- 
  cnstrnt_dta_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::group_by(dta_src) %>% 
  dplyr::summarise(n = sum(n)) %>% 
  dplyr::mutate(l = paste0(dta_src, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = "\n")
# plot base
plt_fshed_cnstrnt_lvl_temp <- ggplot(data=cnstrnt_dta_temp) +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(dta_src)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    # , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
# plot with large labels
plt_fshed_cnstrnt_lvl_temp + 
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.065,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 4 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.065,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.5 # 2.8
    , vjust = 2.3
  ) 

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_dist_sc.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

4.8.3.1 Combine map and fireshed plot

# plot small for inset
inset_temp <- plt_fshed_cnstrnt_lvl_temp +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.105,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3.2 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.105,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.1 # 2.8
    , vjust = 2.3
  ) +
    cowplot::theme_minimal_grid(9) +
    theme(
      legend.position = "none"
      , plot.background = element_rect(color = "black", fill="white", linewidth = 2)
      , axis.text.x = element_blank()
      , axis.title.x = element_text(size=7)
      , strip.text = element_text(color = "black", face = "bold", size = 9)
    )

# combine
cowplot::ggdraw(
  plt_fpas_sc_map[[2]] + 
    theme(
      legend.position = c(.4,.98)
      , legend.direction = "horizontal"
      , legend.text = element_text(size = 9)
      , plot.background = element_rect(fill="white", color = NA)
      , plot.margin = margin(0, 5.1, -0.7, 0.1, "cm")
    )
  ) +
  cowplot::draw_plot(inset_temp, x=.5, y=.56, width=.45, height=.4) +
  cowplot::draw_plot_label(
    label = c("A", "B"),
    x = c(0.02, 0.5),
    y = c(1, 0.95),
    size = 17,
    fontface = "bold",
    color = "gray11"
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_sc_map_inset.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

4.9 Landcover Class Area

4.9.1 Overall Landscape

What is the landcover classification of the area within the landscape boundary? This GEE script was used to calculate the area by NLCD landcover class.

# nlcd class groups
nlcd_reclass_df_temp <- data.frame(
    class_id = c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95) %>% as.character()
    , class_label = c(
      "Water"
      , "Water"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Barren"
      , "Forest"
      , "Forest"
      , "Forest"
      , "Shrubland"
      , "Shrubland"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Planted/Cultivated" # double-check name
      , "Planted/Cultivated"
      , "Wetlands"
      , "Wetlands"
      )
  ) %>% 
  dplyr::left_join(
    FedData::nlcd_colors() %>% dplyr::rename(class_id=ID) %>% dplyr::rename_with(tolower) %>% dplyr::mutate_all(as.character)
    , by = dplyr::join_by(class_id)
  )
# load data by wf priority 
wf_nlcd_area <- readr::read_csv("../data/wildfirepriority_WCS202308/wfpriority_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::left_join(
    wf_landscapes %>% 
      dplyr::select(name,state,area_name)
    , by = dplyr::join_by(name,state)
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )
# get color palette
nlcd_colors <- wf_nlcd_area %>% 
  dplyr::group_by(class_id,color,class_label) %>% 
  dplyr::summarise(sum_area = sum(area_m2,na.rm = T)) %>% 
  dplyr::arrange(class_label,desc(sum_area)) %>% 
  dplyr::group_by(class_label) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::pull(color)
# aggregate and plot
wf_nlcd_area %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    pct_area = sum_area_ha/sum(sum_area_ha)
    # , area_name = area_name %>% forcats::fct_rev()
  ) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "Percent of Overall Landscape Area"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_nlcd_area_ls.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )